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Abstract. We present extensive computational results for the effective temperature, defined by the 
' fluctuation-dissipation relation between the mean square displacement and the average displacement of 

, grains, under the action of a weak, external perturbation, of a sheared, bi-disperse granular packing of 

compressible spheres. We study the dependence of this parameter on the shear rate and volume fractions, 
the type of particle and the observable in the fluctuation-dissipation relation. We find the same tempera- 
ture for different tracer particles in the system. The temperature becomes independent on the shear rate 
for slow enough shear suggesting that it is the effective temperature of the jammed packing. However, we 
also show that the agreement of the effective temperature for different observables is only approximate, for 
very long times, suggesting that this defintion may not capture the full thermodynamics of the system. On 
the other hand, we find good agreement between the dynamical effective temperature and a compactivity 
' calculated assuming that all jammed states are equiprobable. Therefore, this definition of temperature may 

, capture an instance of the ergodic hypothesis for granular materials as proposed by theoretical formalisms 

for jamming. Finally, our simulations indicate that the average shear stress and apparent shear viscosity 
follow the usual relation with the shear rate for complex fluids. Our results show that the application of 
shear induces jamming in packings whose particles interact by tangential forces. 
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1 Introduction 

Granular materials are collections of macroscopic particles 
which interact, essentially, through inelastic collisions, re- 
pulsive contacts, and friction. This idea is simple and can 
be applied to a great variety of systems in nature, such as 
grains of sand, and assemblies of balls. Their importance 
is enormous to some industry segments, for example, the 
food industry [Q. 

Even though granular materials may sound simple, 
their physics is not. For reviews on general properties of 
such systems refer to the work by Jaeger and collaborators 

The complex properties of granular systems make the 
development of a general theory describing their behav- 
ior a difficult problem. A statistical mechanics approach 
was proposed by Edwards [HE] to describe dense, slowly 
driven packings of rigid grains. The strong assumption 
of Edwards is that the microscopic configurations, called 
blocked or jammed, that the system visits during its evo- 
lution are equally probable. By this hypothesis, Edwards 
derived a whole thermodynamics-like formalism for pow- 
ders. 

Makse and Kurchan [S] showed, for system of driven, 
compressible grains, that the quantity: 

where S is the entropy of the powder and E, its energy, 
is the same, within error bars, as an effective temperature 
Tf.f f obtained dynamically, as proposed by Cugliandolo et 
al. ini- The quantity Xe is the compactivity analogous to 
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the compactivity of Edwards T where instead of using the 
energy E it is used the volume V of the powder. This effec- 
tive temperature, Tf,f / , is defined through the fluctuation- 
dissipation relation (FDR) between the mean square dis- 
placement (MSD) and average displacement (AD) of grains 
in z-direction under the action of a weak external force /: 

(W.)-.(0)f)=2T.,,Mi^, (2) 

the brackets indicate an average over time and particles. 
This parameter is also called configurational temperature 
for it controls the structural relaxation of a jammed sys- 
tem. This idea was studied in a variety of systems includ- 
ing glasses PITIISIIMTU] , models of powders QT] and foams 

m 

The purpose of the present paper is to deeper investi- 
gate this result. We study the dependence of the configura- 
tional temperature with the external parameters (volume 
fraction and shear rate). We show that it does not depend 
on the shear rate (for slow enough motion) and the par- 
ticular species of particle, a result that we believe charac- 
terizes the jamming, and provide initial evidence that the 
effective temperature is a proper thermodynamic variable 
to describe the system. We also discuss the dependence of 
Tf.f f on the observable used in the fluctuation-dissipation 
calculation. It has been recently shown that the effective 
temperature of dense granular matter can be measured 
experimentally in a Couette cylinder filled with spherical 
grains and slowly sheared ^31 ■ These results together with 
the work performed by other groups provide indications 
that the thermodynamics approach to dense granular mat- 
ter is possible. 
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We also analyze the role played by the granular tern- 2 Simulation 

perature, as mesured by the average of kinetic energy. The 

We use Molecular Dynamics (also called Distinct Element 

results indicate that it is not a thermodynamical param- 

Mcthod (DEM) [HI) to model our system. It is a bi- 

eter for dense regimes due to its strong dependence on 

disperse assembly of 1000 spheres, half large, half small, 

the external parameters. We also observe the behavior 

with radii ratio of Rs/Rl ~ 0.818, which are initially 

of others macroscopic parameters of the packings, such 

randomly generated in a periodic cubic cell of size L. 

as the average shear stress and apparent shear viscosity. 

The spheres interact via normal and tangential forces 

The first observable follows the usual Herschel-Bulkeley 

|15j . The first one is given by Hertz law, F„. The second 

law. There seems to be a violation of this law for states 

force is Mindlin's no-slip solution for the contact. Ft- The 

in low densities where there is a frictional force between 

tangential force can increase up to the Coulomb threshold 

grains. Besides, the denser states develop a yield stress 

of friction Ft = jiFn (where fi is the friction coefficient, 

at low shear rates. These dense states are in the regime 

typically /i — 0.3). If Ft increases beyond /i-F,;, the grains 

where the effective temperature is the relevant parameter. 

slide and the contact is broken. We also simulated packings 

This implies the same correspondence already observed 

in which particles interact only through normal forces, and 

for glass-forming liquids 0^- The apparent shear viscos- 

we refer to these cases no friction, or frictionless, systems. 

ity displays the shear-thinning behavior of complex fluids. 

We notice that the frictionless case would model a system 

suggesting that powders may be treated as such fluids. 

of compressed emulsions Gravity is absent in all 

our simulations. 

During the generation phase, no overlap between par- 
ticles is allowed. We commonly have a initial volume frac- 
The organization of this work is as follows. In Section tion, 0, of the order of 0.25 in all our simulations. Since 
II we present details of the procedures employed in the our goal is to study dense systems, we should compress 
simulations. In Section III, the results for the diffusion the initial packing in order to attain a higher volume frac- 
coefHcients, mobilities, deflned by the MSD and AD, re- tion. It is known ^Tj that there is a jamming transition 
spectively, and effective temperatures are presented. Also at a critical volume fraction 0c, corresponding to the ran- 
in this section, we discuss the role of the granular temper- dom close packing (RCP). Above RCP there exists a me- 
ature. The results of the average macroscopic quantities chanically stable solid-like state which is characterized by 
measured, shear stress and apparent shear viscosity, are a non-vanishing internal stress and coordination number, 
shown in Section IV. Finally the conclusions are drawn in Therefore, we choose three volume fractions, cf) = 0.6428, 
Section V. 0.6623, 0.7251, above RCP to perform our studies. 
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We study the jamming transition by performing sim- 
ulations also below RCP, i.e., in the fluid-like state. Such 
states can be prepared by expanding a previously equili- 
brated system above RCP to decrease the volume fraction 
up to the point where the pressure and the coordination 
number vanish, as reported in |I7] . We prepare three fluid- 
like states with volume fractions given by cf) — 0.5999, 
0.6104, 0.6318. 

The simulations involve the application of shear in the 
x-y plane, at constant volume (see Fig. with x the 
direction of the flow and y the velocity gradient. We follow 
the trajectory particles in the z-vorticity direction. We use 
a modified version of the usual Lees-Edwards boundary 
conditions |18| which imposes a linear velocity profile in 
the shear plane. Normal periodic boundary conditions are 
enforced in the z-direction. Segregation was not observed 
in the time scale of our runs. 

For each volume fraction, 9 different shear rates were 
considered, spanning a total of 5 decades of magnitude: 
7 = 10"^ 5 X 10"^ 10"^ 5 x W-^, m-^, 5 x lO^-*, IQ-^, 5 X 
10"'^, 10~^ s^^. In our simulations time is given in seconds, 
length in meters, and mass in kilograms. 

In what follows, we will calculate all the observables 
above and below RCP and for frictional and frictionless 
particles so that we will investigate the jamming tran- 
sition at the RCP limit. The results shown in the next 
sections will be only for large particles, unless otherwise 
noticed. The results will be presented as functions of the 
shear rate, one curve for each volume fraction where appli- 
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a) Siinulations 




Fig. 1. Detail of the simulations of grains of 100 microns in- 
teracting via Hertz-Mindlin contact forces. A slow shear flow, 
indicated by the horizontal arrows, is applied to the jammed 
system. We follow the tracer particle trajectories to obtain the 
diffusivity. An external force F is then applied to the tracers 
in response to which we measure the particle mobility. These 
dynamical measurements yield an "effective temperature" ob- 
tained from an Einstein relation. 

cable, in both frictional and frictionless cases. Moreover, 
since sheared systems enter a stationary state after an ini- 
tial transient time jl9ll2()j and does not display aging, it is 
sufficient to calculate all the following quantities as sim- 
ple time averages. We calculate correlation and response 
functions with sets of overlapped data in the same run. 
We call this procedure window average. In practice it is 
the same as taking several time series in different runs and 
averaging over them. It saves run time and storage space. 

3 Effective temperature 

We calculate the effective temperature, T^ff, of the gran- 
ular system by measuring it through eq. (jSJ, as in 0. We 
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investigate its dependence on the shear rate, and volume guides to the eyes, unless otherwise noticed. Figure llfa) 
fraction. Next, we show that Te// does not depend on the shows the MSDs of the particles versus time from where 
particular species on which it is being measured. Finally, the diffusion constant is extracted. 

we measure Te// by a different FDT commonly used in In Fig. Oand^] we show the results for for fric- 
the glass literature, the self-intermediate scattering func- tionless and frictional cases, respectively. First, we notice 
tion (ISF) and averaging over the jammed configurations, that the high volume fraction curves go to zero in a similar 

fashion, suggesting that our system is jammed |21| . Also, 
as the volume fraction decreases, particles have a higher 
diffusion for smaller shear rates. This is a consequence of 
the larger space available for particles to move under the 
influence of enduring contacts. As we go up in shear rate, 
this difference tends to disappear, the diffusion coefficient 
becomes independent on the density, due to the larger ef- 
fect of inter-particle collisions. 



time (10 cycles) 




Fig. 2. (a) MSD and (b) AD for large (lower) and small (up- 
per) particles in units of large particle diameter, di . 



3.1 Diffusion, mobility and FDT 

The calculation of the diffusion coefficient in the z-direction, 
Dz, requires the evaluation of the MSD in this direction 
as a function of time: 

N 

N 




0.7251 
0.6623 
0.6428 
^0.6318 
^ 0.6104 
0.5999 



/ 1 ^ 



(3) 
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Fig. 3. Diffusion coefficient versus shear rate, no friction. 



In the case of friction, particles still diffuse at lower 



where N is the number of particles considered, Zi(t) is the densities but, somewhat surprisingly, the diffusion coeffi- 

z-coordinate of the i-th particle at time t. The diffusion cient becomes independent on the density even for (p = 

coefficient was calculated separately from the mobility (see 0.6318, a fluid-like state. This is a first indication that the 

below). In all plots, lines connecting the points are only presence of friction between particles can induce the jam- 
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ming transition for systems that are below RCP. This is 
usually called a shear-induced jamming transition, see |22l 
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Fig. 4. Diffusion coefficient versus shear rate, friction. 

The mobility, x^, is the response function associated 
with the motion of particles under an applied external 
force. To calculate Xz? ^ small constant force is applied 
to each particle in order to induce a displacement in the 
z-direction. The mobility is taken as the slope of the curve 
of the average z-displacement, normalized by the applied 
external perturbation, versus time: 



(4) 



where / is the applied force, and = ±1 is a "charge" 
set to each particle i. This is a standard procedure em- 
ployed in numerical calculations to improve averages |3: • 
An average over a distribution of these charges yields the 
desired property. 




time (10 cycles) 



Fig. 5. Average z-displacement (given in units of large particle 
diameter) for force amplitudes of (from top to bottom) 0.06, 
0.05, 0.04, 0.03, 0.02, 0.01, 0.009, 0.008, 0.007, 0.006, 0.005, 
0.004, 0.003 N. 

Before calculating the mobility, we make short runs 
in each state considered for, at least, 10 forces in order 
to measure {z{t) — z(0)) //. The linear response regime is 
the range of forces where all z-displacement curves, nor- 
malized by the applied force, should overlap within nu- 
merical error. An example of this procedure is shown in 
Fig- El Here, we show several displacement curves, for dif- 
ferent applied forces, in the case of large, frictionless par- 
ticles at volume fraction oi cj) — 0.6428 and sheared at 
7 = 5 X 10^® s^^. We can see clearly that for small forces, 
all [z{t) — z(0)]// curves (without averaging) are within 
some well defined range, while at higher forces, the curves 
starts to deviate from the linear regime. After this regime 
is identified, we selected at least 5 values of the force to 
perform the calculation in longer runs. 

Figs. El and [71 show the results for Xz for frictionless 
and frictional cases. Similarly to the case of diffusion, the 
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Fig. 6. Mobility vs. shear rate, no friction. 

mobility of the grains increases with decreasing volume 
fraction. Also, the presence of friction between grains di- 
minishes the ability of particles to move, hence yielding 
a smaller mobility for these packings. The — 0.6318 
friction packing has similar dependence on 7 as its denser 
counterparts, confirming the suspected shear-induced jam- 
ming transition indicated in the diffusivity plot of fig. 0] 

A clear trend observed in the plots for and Xz is 
that, at lower volume fractions and shear rate, both 
and Xz have weaker dependence on shear rate, clearly ob- 
served in the no friction plots (and particularly, in the 
(/) = 0.6428 and = 0.6318 cases). This fact indicates 
that these states are not jammed, since the ratio between 
Dz and Xz, and consequently T^ff, would not be indepen- 
dent on 7. A similar trend is displayed in the friction plot, 
but only to the two states at (p = 0.5999 and (p — 0.6104, 
and for shear rates above 7 = 10~^ s~^. 

At higher volume fraction, we have a different picture. 
The dependence on the shear rate is approximately the 
same for all curves in all plots. This already suggests the 



regime where Tef / would be the relevant thermodynamical 
parameter. In these figures we put a dashed line with unity 
slope in order to indicate the expected behavior of both 
Dz and Xz with the shear rate. We have close agreement 
in our results. 
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Fig. 7. Mobility vs. shear rate, friction. 

In Fig. we plot the MSD in the z-dircction while 
in Fig. 133, the AD in z for large and small particles after 
the external force is applied. We find that the small tracer 
gives rise to a significantly larger diffusion constant than 
the large tracer. The same behavior is obtained for the 
mobility as can be seen in Fig. This case was the (p — 
0.7251, frictionless one sheared at 7 = 5 x 10^^ s^^. The 
perturbation was / — 0.09. All lengths are in units of large 
particle diameter, di. In fig. we measured Dzld\ = 
2.295 X lO"'^ s^\ large, and Dz/dl = 3.995 x lO"'^ s-\ 
small particles. Figure gave Xz/di — 2.232 x 10^® 
s/(kg-m), large, and Xz/di — 3.741 x 10^^ s/(kg-m), small 
particles. The diffusion constant and the mobility are in- 
versely related to the grain sizes. For a Stokes fluid we 
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would expect that the ratios of the mobihties and diffu- 
sivities satisfy Dg/Di — di/dg, and Xs/xi — di/ds where 
ds and di are the diameter of the small and large particles 
respectively. 

In our simulations, the ratio di/dg is equal to 1.22. We 
have, for example, the diffusion coefficients for the — 
0.7251 friction case, 7 = 10"^ g-i^ ^re A = 9.317 x 10"" 
and Ds = 1.149 x 10"^^ their ratio yields 1.233, thus 
in good agreement with the Stokes result. Moreover, the 
mobilities for the state at (/> = 0.6104 volume fraction, 
sheared at 7 = 5 x 10^^ s^^ without friction have the 
values xi = 6.818 x 10^* and Xs = 8.664 x 10~^, whose 
ratio gives 1.27. 




10 10 

Y(l/s) 



Fig. 8. Effective temperature vs. shear rate, no friction. 

Now we pass on to the results of effective tempera- 
ture. They are given in Figs. |S1 and |51 We can see that 
at higher volume fractions this effective temperature is 
approximately independent on the shear rate as 7 — > 0. 
This confirms what was stated before about the range of 
jammed states. Also, it is an important result for it shows 



that this is indeed an intrinsic property of the granular 
system, since it does not depend on the driving strength. 
The jammed state can be characterized by the indepen- 
dence of the effective temperature on the shear rate since 
this behavior is not observed in the other cases studied. 




y (i/s) 

Fig. 9. Effective temperature versus shear rate, friction. 

Two other conclusions may be drawn from these graphs. 
One is already expected from the previous plots: the state 
— 0.6318 for the frictional case is indeed jammed, al- 
though it is below RCP. This confirms that the application 
of shear can induce a jamming transition in a friction pack- 
ing. The second conclusion is that the state (/) = 0.6428 
for frictionless particles is not jammed, even though it is 
above RCP. This indicates that frictionless particles are 
unstable under shear, however small. 

Next, we test the zeroth law of thermodynamics for 
Tg//. In other words, if the effective temperature has a 
physical thermodynamic meaning, it should be the same 
for different tracers in the medium. We test this by re- 
peating the calculations of the MSD and the AD for both 
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types of particles in frictionless states (/) = 0.6428 and 
(t> — 0.6623. Figure [TUI shows the parametric plot of the 
two quantities at the cases considered, sheared at 7 = 
5 X 10^^ s^^. The agreement is clearly excellent. Although 
the curves for (f) = 0.6428 are in close agreement, the re- 
sults for Teff rules out the possibility that this state is 
jammed. This further suggests that Te// is a good ther- 
modynamic variable for the system. 




— 0.6428, Large, no fric 
-- 0.6428, Small, no fric 

— 0.6823, Large, no fric 
■ ■ ■ 0.6823, Small, no fric 
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Fig. 10. Parametric plot of the MSD by AD. The legend dis- 
plays the cases studied. All lengths are given in units of large 
particle diameters. 



3.2 Self-Intermediate scattering function FDT 

Another important characteristic of the effective temper- 
ature is that it should be independent of the particu- 
lar observable (or fluctuation-dissipation relation) used in 
its measurement. In other words, we would like to know 
whether by measuring Tf,ff using different observables 
than those in |(5J) we obtain the same result. This inde- 



pendence on the observable was showed for supercooled 
liquids ^U] and foams \V2\ - 

Here we use a quantity which was extensively studied 
in the glass literature, namely the self-intermediate scat- 
tering function (ISF) |51l^lOII25ll26ll27j . defined as: 



CkAi) = (^I]exp{ifc, [z,{t)-z,{m 



(5) 



1=1 I 



(6) 



where is the wave vector in the z-direction. This func- 
tion is the Fourier transform of the self part of the Van 
Hove correlation function at time t |211. We chose the wave 
vector kz as the peak of the static structure factor, which, 
in our case, is k^ — 71.8L~^. 

In order to measure a FDR 10 we apply a small force 
in each particle i of the form Fi — —fkz^i sm{kzZi)^ where 
/ is the amplitude of the force. After the perturbation is 
applied we measure the response function, defined as: 

N \ 

f \N 

The amplitude / is the quantity which controls the inten- 
sity of the force. Therefore we set it to be small enough 
so that we measure Rk^ in the linear response regime. 
As in the mobility calculation, we search for the linear 
response regime by measuring the response function for 
several forces before performing the full calculation. Also, 
we did not employ the window average in the calculation 
of lO, instead we averaged over different measurements, 
20, for each force. 

The effective temperature is defined as the ratio be- 
tween the real parts of ^ and © : 



Re[RkM' 



(7) 



10 
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3.3 The role of kinetic energy 
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Fig. 11. Tfie effective temperature (straigtit line) calculated 
from ^ overlaps with fluctuation-dissipation plots from Q 
in the very long time limit. The amplitudes of the applied 
perturbations are (from top to bottom) 5 x 10~''', 6 x 10~^, 
7 X 10~^ and 8 x 10"^ N. This case is for frictionless particles, 
volume fraction of <^ = 0.6428 and shear rate of -y = 5 x 10~^ 

We measured this effective temperature for the (p = 
0.6428 frictionless state sheared at 7 = 5 x 10^^ s^^. In 
order to test the equivalence of both temperatures, we 
plotted all results for ^ parametrically and a straight 
line with slope —Teff~^, with its value given by the cor- 
responding result from ^ . This is shown in Fig. 1111 for 
large particles. We see that the straight line has slope in 
the order of magnitude, only for very long times, as the 
one obtained from the ISF. Thus, it possibly indicates 
that Teff do not capture all the statistical mechanics of 
the jammed state. 



The kinetic energy is defined as the sum of the velocity 
fluctuations of all particles: 

N 

i=l 

where uii is the i-th particle's mass and , is the velocity 
vector. 

As in thermal systems, we can define a parameter, 
called granular kinetic temperature, Tk, which is the aver- 
age of JHl divided by 3A^. Since the average velocity is zero 
in the z direction, since there is no bias in this direction, 
we take Tr- as the following expression: 

1 ^ 

i=l 

This is an estimate of the temperature of the fast modes 
of the system [Hlll9j (since a granular system is athermal 
by definition). In thermal systems, eq. © yields the bath 
temperature. It is expected, in the dense limit, that Tk 
does not share the same characteristics of T^f / . 

In Fig. ^1 we show the kinetic temperature for both 
large and small particles for the highest and lowest vol- 
ume fractions studied in the frictionless case. It is directly 
seen that such plots evidence the difference of Tk for both 
species of particles. We also notice that the kinetic tem- 
perature has two different regimes. One at high shear rate, 
where it is almost density independent. This is the fluid 
regime where the kinetic energy is dissipated at a smaller 
pace than in the dense case, this dissipation made mainly 
by collisions between grains. At low shear rates, the ef- 
fect of volume fraction shows up. The denser the system. 
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the smaller the kinetic temperature, as expected, since at 
high volume fractions there is a large number of contacts 
in the packing, which rapidly dissipate kinetic energy. In 
this case, the kinetic temperature depends on the shear 
rate again as a power law, but with a density dependent 
exponent. 




Fig. 12. Average kinetic temperature vs. shear rate, no fric- 
tion, for large and small particles. The legend displays the cases 
showed. The dashed lines are fits to the curves and the numbers 
are the slopes of these lines (explained in the text). 

Figure [TSl shows the kinetic temperature for large and 
small particles, in the lowest and highest volume fractions 
considered here for frictional packings. The same conclu- 
sions drawn from plot 1 121 can be taken here. This parame- 
ter is indeed different for both particle species, confirming 
that Tk should not be used as a thermodynamical pa- 
rameter for dense systems. As in the previous graph, Tk 
has two different regimes. One for high shear rate, where 
it scale with 7 as a power law with exponent 1.26. The 
second one, where the jammed states are, this exponent is 



density dependent, the lowest being 0.61 and the highest 
1.01. 
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Fig. 13. Average kinetic temperature versus shear rate, fric- 
tion for large and small particles. The legend indicates the 
cases showed. The dashed lines are fits to the curves and the 
numbers are the slopes of the lines (explained in the text). 

All dashed lines and their respectively slopes are fits to 
the curve of Tk for the total number of particles, without 
any distinguishing from large or small. It should be no- 
ticed that the results of 28 indicate that the kinetic tem- 
perature is linearly dependent on the driving strength, the 
constant of proportionality being a function of the density 
and elasticity of the components of the mixture. Here we 
obtained a power law dependence between the two quan- 
tities. It should be noticed, however, that the previous 
results are for a granular bed, while ours is a 3D system. 
This difference may induce different dependencies in these 
parameters. 

A comparison between FigureslT^andlTSIwithlSlandl^ 

respectively, implies that the effective temperature is al- 
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ways larger than the granular one in the jammed regime. 
We can also see that the behavior of Tg// with volume 
fraction is different from the behavior of Tk with 0. The 
last one decreases with increasing while the former in- 
creases. Also, one can clearly see that the ratio of the val- 
ues of Tk for large and small particles is also constant, in- 
dependent on the shear rate. This is the same observation 
that was made at the experiments of j28ll29j . Their result, 
which was obtained for a system in a fluid-like state, also 
applies to a dense one, at least to the component T^z of 
the temperature tensor (as defined in 30 ). We did not 
studied the other components of this tensor. Finally, we 
see that the kinetic temperature is always larger for large 
grains, a result also observed in [2511^ . 



The growth of the exponent of Tk in 7 with volume 
fraction can be seen in terms of the interaction between 
grains. Shearing a packing puts energy into it and the 
interactions dissipate this energy. As we lower the shear 
rate, enduring contacts build up as the main interaction 
mechanism. Since at lower volume fractions there are less 
contacts, diminishing 7 does not cause a sharp drop on 
the average energy. For denser cases, the same reduction 
in shear rate causes the average energy to greatly lower 
its value due to the presence of more contacts. Therefore, 
Tk drops faster with increasing <j). Moreover, exponents are 
larger for friction packings due to the extra dissipation 
provided by tangential forces. 
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4 Calculation of the temperature via 
configurational averages 

So far we have shown the existance of an effective temper- 
ature for jammed granular materials and argue that this 
could be a thermodynamical variable of the system be- 
cause two different subsystems characterized by particles 
of different sizes have the same value of this parameter, as 
shown in fig. ^] This is a kind of zero-th law for the sys- 
tem. This property is not shared by, for instance, the gran- 
ular temperature, and therefore it could not be thought 
of as a physical variable describing the thermodynamics 
of grains. 

We also showed that the effective temperature is within 
the range of the temperature measured by other observ- 
ables but only for the very long time behavior of the sys- 
tem, as shown in Fig. 1111 Thus our results suggest that 
this temperature may not capture all the thermodynamics 
of jamming. 

Another test of the physical significance of Tg/ / is to 
compare its value with a temperature obtained directly by 
an average over the jammed configurations. 

If it was true that a thermodynamic framework could 
describe the behavior of jammed systems, it stands to 
reason that the compactivity Xe of the granular packing 
given by eq. should be measured from a dynamical ex- 
periment involving the exploration of the jammed states. 
If the system is ergodic and amenable to a statistical me- 
chanical approach, then the temperature defined by the 
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configurational average should be the same as a dynami- 4.1 Exploring the jammed configurations via a flat 
cal measurement, as done for instance by the FDT. average. Test of ergodicity: Tefr — Xe 



To test this assumption is to test the ergodic hypoth- 
esis for jammed matter. Indeed, the statistical framework 
proposed by Edwards 3 is based on the fact that the 
jammed configurations arc all cquiprobable. This assump- 
tion is analogous to the basis under which the micro- 
canonical ensemble of equilibrium statistical mechanics is 
built. However, this assumption has been very controver- 
sial for granular matter since grains are path dependent, 
frictional and non-conservative. In the absence of a Hamil- 
tonian to describe the dynamics, there is no Liouvillc's 
theorem on which we can construct the statistical me- 
chanics of grains from first principles. Thus, the idea of 
ergodicity for granular materials has received much oppo- 
sition. 



The independent method to study the configurational space 
of grains allows us to investigate the statistical proper- 
ties of the jammed states available at a given energy and 
volume. In turn, it is investigated whether it is possible 
to relate the dynamical temperature obtained above via 
a diffusion-mobility protocol to the configurational com- 
pactivity based on the flat average over jammed states. 

In order to calculate Xe and compare with the ob- 
tained Tcff, it is needed to sample the jammed configu- 
rations at a given energy and volume in an equiprobable 
way. In order to do this, the jammed configurations are 
sampled with the following probability distribution: 

P. ^ exp[-i?7T* - E^_JT,^,] (10) 



The validity of this statement was analyzed through Here the deformation energy E corresponds to the Hertzian 
computer simulations by Makse and Kurchan 5 . Our present energy of deformation of the grains. The extra term added 
results are an extension of these results. Here we reproduce in Eq. (jTUl) allows the flat samphng of the jammed states 
the discussion of this paper to make an explicit conncc- to be perfomcd. The jammed energy is such that it van- 
tion with previous work. We will calculate the temperature ishes at the jammed configurations: 



-Ejammcd OC ^ , (II) 



from the entropy of the packing with a mathematical con- 
struction which assumes the jammed states to be equally 
probable. The temperature calculated from the entropy where Fi is the total force exerted on a particle by its 



of the packing is compared to the effective temperature neighbors. The main property of this energy is that it 

obtained dynamically. This is achieved via the procedure should be zero at jamming. 

applied above which has for an aim to probe each static Two "bath" temperatures are introduced (these tem- 

configuration by allowing the system to evolve at a very peratures are ~ lO^** times the room temperature) which 

slow shear rate. This is done through the FDR which does allow the exploration of the configuration space and the 

not assume the flat average. calculation of the entropy of the packing assuming a flat 
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Fig. 14. Annealing procedure to calculate Xe at different 
elastic compressional energies. We plot the elastic energy vs 
Taux during annealing together with the distribution of elas- 
tic energies obtained during shear (dashed curve, mean value 
(E) = 8.4 X 10""'). At the end of the annealing we find that 
the energy is very similar to the average energy for the sheared 
system only when we set T* — Xe- 

average over the jammed configurations. Equilibrium MD 
simulations are employed with two auxiliary "bath" tem- 
peratures (r* , Taux) , corresponding to the distribution pOfl . 
Annealing Taux to zero selects the jammed configurations 
(-Ejammod = 0), whilc T* fixcs the energy E. 

In practice, the equilibrium MD simulations are carried 
out with a modified potential energy: 



U = I^E + E- 



'jammed; 



(12) 



and calculate the force on each particle from F ~ — VC/. 

The auxiliary temperature Taux is controlled by a ther- 
mostat which adjusts the velocities of the particles to a ki- 
netic energy determined by Taux- We start by equilibrating 
the system at high temperatures (Taux and T* oo) and 



anneal slowly the value Taux to zero and tune T* so as 
to reach the value of E that corresponds to the average 
deformation energy obtained during shear. 
The partition function is 



Z = ^cxp[-i?7T* - E^_JT,^^], 



(13) 



from where the compactivity reads: 



rri* _ dE Taux ^ „ 

^ -dS ^ 



(14) 
^ 0), 



Thus at the end of the annealing process (Taux 
T*{E) = Xe{E), since in this limit the sampled configu- 
rations have a vanishing number of moving particles at a 
given E. 

Since a calculation of the force from a potential en- 
ergy, eq. (|12|l . is required, only conservative systems can 
be studied with the method. Therefore it is not possible to 
test these ideas in a granular system with friction for it is 
a path-dependent system. Thus we focus our calculations 
on a system of frictionless particles. 

A system of frictionless viscoelastic spherical particles 
could be thought of as a model of compressed emulsions 
|31ll32j . see also [33] for details. Even though they can be 
modeled in this way, an important difference arises in the 
inter-droplet forces, which are not given in terms of the 
bulk elasticity, as they are in the Hertz theory. Instead, 
forces are given by the principles of interfacial mechanics 
|34| . For small deformations with respect to the droplet 
surface area, the energy of the applied stress is presumed 
to be stored in the deformation of the surface. The sim- 
plest approximation considers an energy of deformation 
which is quadratic in the area of deformation 34 , analo- 
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gous to a harmonic oscillator potential which describes In the calculation of the conipactivity for the same 
a spring satisfying Hooke's law. More elaborated mod- system of 200 frictionlcss particles at the same pressure 
els have been proposed, and in a recent study the force and volume fraction, there is no dynamics so the viscous 
law was calculated experimentally using confocal micro- forces of dissipation are completely disregarded. At the 
scopic images of emulsions |16| . Here, we will avoid issues end of the annealing protocol the compactivity at a given 
of path-dependency introduced by the interparticle tan- deformation energy can be obtained as illustrated in Fig. 
gential forces, by excluding these forces from the calcu- El 
lation {Ft = 0). Because there are no transverse forces, 

In the begining, the system is equilibrated for 40 x 10^ 

the grains slip without resistance and this procedure es- 

iterations at high auxiliary temperatures, such as at point 

sentially mimics the path to jammed states for the com- 

A in Fig.ini (T* = 3.4 x lO^^, Taux = 3 x IQ-^). Different 

pressed emulsion system. 

annealing paths are taken for T*, one decreases its value 

„, ,TT^ 1 1 n , . • 1 • 1,1 to T* = 2.9 X 10^^ (open squares, path la) and the other 

Ihus, the MD model oi granular materials is adapted 

., , ^ , , • , , anneals it to T* = 3.4 X 10^^ (filled circles, path 1). Path 

to describe the system of compressed emulsions by only 

... , , , r , ;• 1 1 ■ ■; 1 la is followed until the value of T* reaches 2.9 X 10^^. Path 

excluding the transversal forces (tangential elasticity and 

„,,„.. ^ „i • 1- • 1 1 • 1 1 1 1 continues until both temperatures are annealed to the 

Coulomb irictionj . i he continuous liquid phase is modeled 

. .; . 1 r ■ in point where this path is split again: one branch anneals T* 

m its simplest form, as a viscous drag lorce acting on every 

, , . 1 , , . mi 1 to 2.9 X 10^^ (open diamonds, path lb), ioining the line of 

droplet, proportional to its velocity, ihe systems under 

, • 11 1 1 in path la, while the other anneals to T* = 3.4 x 10^"*. From 

investigation have exponentially large (m the number oi 

. , , 1 r 1 1 • 1 11 this point on, T* is kept constant in paths la and lb, while 

particles) number of stable states jammed at zero bath 

„, . n 1 1- 1 • annealing Taux to zero. Path 1 continues until it reaches 

temperature, thus, a meanmgiul annealing procedure is 

, -1 1 r 11 ■ -11 point B, with Taux = 3 x 10^^*^), where it is split for the 

only possible tor very small system sizes with the present 

, mi r 1 nil l^-st tlmc 111 thrce different branches: one keeps T* constant 

computational power, iheretore, the system ot study has 

,-11 at 3.4 X lO^"' (filled right triangles, path 3) , the second 

200 particles only. ^ or-/, 

anneals T* to the value 2.9 x 10^''^, the last one anneals 

The calculations for the effective temperature for this T* to 3.4 x 10^^ (filled triangles, path 2). When T* = 

system have been done in 5, . The effective temperature Tcff (Paths 1, la and lb), the final elastic compressional 

is Teff — 2.8 X 10^^, the mean pressure P = 10 MPa, energy value when Taux falls inside the distribution of 

the volume fraction (f> = 0.66 and the mean compressional energies obtained, and it is very close to the mean value 

energy (E) — 8.4 x 10^"* J. All these values have been of the elastic energy during shear (E). The remarkably 

obtained during the shear of the material. result is that this implies that the compactivity and the 
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effective temperature obtained dynamically are found to 
coincide to within computational error, 

Xe « Teff. (15) 

For other values of T* ^ Tcs the final E falls out of 
the distribution obtained during shear (Paths 2 and 3). 

This provides evidence for the validity of the effective 
temperature as a dynamical estimate of the compactivity, 
and more importantly, justifies the use of the statistical 
measurements presented in characterizing the macroscopic 
properties of the jammed system. Since T^ff does not as- 
sume the flatness of the jammed states because it is ob- 
tained dynamically, these results indicate that, at least in 
the present system, the ergodic hypothesis for granular 
matter is a plausible approximation to the dynamics of 
the system. 

The conclusion is that the jammed configurations ex- 
plored during shear are sampled in an equiprobable way as 
required by the ergodic principle. Moreover, the dynamical 
measurement of compactivity renders the thermodynamic 
approach amenable to experimental investigations |13j . 

It should be noted though that this is only an approx- 
imation, and it is valid only when the granular material is 
pretreated so that it is in a reversible state. In computer 
simulations, the reversible state is achieved by simulating 
packings without friction as explained before. In the case 
of friction, it is still possible to generate packings which are 
reversible and amenable to a statistical description. This 
problem has been treated before in |16II35| . Experimen- 
tally, it has been also shown that granular matter should 
show reversibility if it is properly prepared, even when 
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system has friction and path-dependent forces. The sem- 
inal experiments of the Chicago group of column tapping 
|36II37| clearly indicate the existance of these reversible 
jammed states in granular materials. Moreover these ex- 
periments have been reproduced by other groups as well 
|38j and under other conditions such as oscillatory pres- 
sure [35] • These reversible states are the ones for which the 
statistical mechanics is possible. Notice that this assump- 
tion preclude other simple systems such as a pile of grains. 
A sandpile, for instance, is not in a reversible state since 
any perturbation will deform the pile. Moreover, there is 
a whole branch of irreversible states that have been ob- 
tained in the Chicago experiments and in the experiments 
of Brujic et al. jSH], before the reversibility is achieved. 

The present thermodynamics is not applied to such 
systems. Other theories, such as the fragile matter theory 
of Gates et al. 021 might be able to capture this irre- 
versible states of granular matter. 

5 System properties 

In this section, we investigate how other properties of our 
system behave as a function of the shear rate and density. 
The properties observed were the average potential energy, 
pressure, and shear stress during shear (only shear stress 
from now on) and apparent shear viscosity. We want to 
know how jamming affects these properties. We only show 
results for the shear stress and the viscosity, since the 
potential energy and pressure behave qualitatively as the 
shear stress. All the properties are calculated when the 
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system enters the steady state. Therefore, they are given the yield stress cto in H16|l . It is not surprising to observe a 



as simple time averages. 



5.1 Shear stress 

The shear stress is defined as the a;?/-component of the 
stress tensor cTy. In fact, we measured the shear stress 
which arises from compression and do not include the con- 
tribution from momentum exchange between grains. In 
complex fluids, the average shear stress is usually fitted to 
a power law, the Herschel-Bulkeley equation |40| : 



{<7xy) = To + A7", 



(16) 



where ctq is the yield stress and A is a constant [TUIl41ll42| . 
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Fig. 15. Average shear stress vs. shear rate, no friction. The 
plateaus are the evidence of a jamming transition. The dashed 
lines are explained in the text. The numbers are their slopes. 

In Fig. ^1 is shown the results for aij for the friction- 
less case. It is observed a plateau at higher volume frac- 
tions and small shear rates. The plateaus correspond to 



non- vanishing yield stress since equilibration of such pack- 
ings implies that. What is interesting is the absence of an 
yield stress for the = 0.6428 state, initially in the solid- 
like state. This confirms that this system is unstable under 
shear and it is not truly jammed. Also, we let the system 



relax after shearing and observe that, indeed, {axy) 



0. 



The existence of a yield stress can be considered as an 
instance of the jamming transition 
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Fig. 16. Average shear stress versus shear sate, friction. As 
in fig. 1151 the plateaus correspond to jammed states. Dashed 
lines, showed with the correspondent slopes, are explained in 
the text. 

As we go up in shear rate, the stress increases as a 
power law whose exponent depends on the volume frac- 
tion (in the cases of the thermal systems, this exponent de- 
pends on the bath temperature |41II42| '). The upper dashed 
line in E| is a fit to the last 3 points in the flow curve 
for the (f) — 0.7251 case. The lower curves do not display 
a yield stress, but they match the opposite limit of the 
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phenomenology of the shear stress, the power-law-fluid 
regime |19II41II42| . The lower dashed line is a fit to the 
flow curve of the = 0.5999 state, showing that this state 
corresponds to a Newtonian fluid. The two middle curves 
match the behavior of power-law fluid, with the exponent 
n of H16|l density dependent. 

In fig.llfilwe have the frictional results for {(Txy)- We no- 
tice that there are four states that develop an yield stress, 
the same four that display a plateau in the effective tem- 
perature results. As in the frictionless case, at high shear 
rate the main contribution to the shear stress is the sec- 
ond term on right hand side of Hlt)|) . Again the exponent n 
depends on volume fraction. Its lowest value comes from 
the slope of the 3 last points in the (a^y) curve for the 
(/) = 0.7251. On the other hand, it seems that, at lower vol- 
ume fractions, the systems displays two different regimes, 
without an yield stress, in conflict with what is expected 
from the Herschel-Bulkeley equation. These two regimes 
are exemplified by the two lower dashed lines which arc 
fits to the first 4 points and to the last 5 points in the 
(/) = 0.5999 flow curve. We would expect this curve to 
have a unity slope and this to be the Newtonian regime 
This fact deserves more investigation, since friction is in- 
herent to grains, and its presence could alter significantly 
the behavior of the shear stress at low shear rate. The im- 
portant conclusion is that at high shear rate the exponent 
n is still density dependent and increases with smaller den- 
sity, eventually reaching 1. Moreover, we also expect that 
at higher shear rates, the average shear stress should scale 
as {(Txij) ^ 7^, the well-known Bagnold scaling |45j . 
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By looking at the results for T^ff, the plots for shear 
stress confirm the two conclusion drawn before: instability 
of frictionless systems at densities up to (p — 0.64, at least 
(this is evidence of the fragility concept proposed to de- 
scribe granular matter by Gates et al. and the shear- 
induced jamming transition. We can see that the plateaus 
in the shear stress (the yield stress) , in the frictionless and 
frictional cases, correspond to the plateaus in the plots for 
Teff. Therefore, these states are the jammed states of our 
system. In this way, we can say that the appearance of a 
well-defined effective temperature is a consequence of the 
development of an yield stress |10II19II43II44| . 

The packing instability can be explained by the fact 
that when shear is stopped, the grains become less com- 
pressed by pushing one another in order to attain a mini- 
mum of potential energy. While they push each other, they 
also slide related to their neighbors, breaking the contacts. 
When all, or most, of these contacts are broken all quanti- 
ties that depend on the overlap between grains (potential 
energy, pressure, shear stress), should drop to zero, since 
the amount of compression is negligible. This relaxation 
does not happen at higher (p because there is no available 
space for the grains to fully relax. 

In the friction packings, the same mechanism of re- 
laxation is present. This difference in the flow curves is a 
consequence of the fact that the remaining normal forces, 
responsible for separating the particles, are not strong 
enough to break the contacts, which are supported by the 
tangential forces. Therefore, we say that our system un- 
dergoes a shear- induced jamming transition in the case 
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4> = 0.6318. The shear induces contacts between particles 
and these shear-induced contacts survive after the shear 
is stopped, providing the system an yield stress. 



5.2 Apparent shear viscosity 

The last of the macroscopic observables probed is the ap- 
parent shear viscosity of the system, 77, defined as: 
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Fig. 17. T) normalized by Tk^^^, frictionless case. The slope of 
the dashed line indicate the regime we can find the system. 

shear rate and the kinetic temperature are coupled. 
Since this input in energy leads to higher T/^, more col- 
lision take place and, consequently, the system reaches a 

higher stress. We should, then, plot the viscosity normal- 

1/2 

ized by T^' in order to disregard this effect on 77. 



Using eq. (|16|l . we can rewrite this expression as 

r] = Aj-'^'-'^l (18) 

From this equation we can see that, according to the value 
of n, the system may display three different regimes: shear- 
thinning {a{4>) > 0), Newtonian {a{(j)) — 0) or shear- 
thickening (a((/)) < 0). Measurements of viscosity have 
been performed in a variety of complex fluids, for example 
glasses |20II19| . foams |^TII^ . and colloids (47II23II48II49| . 
In some of these systems, a shear-thinning behavior, with 
exponents between 0.6 and 1.0, was observed |19II20II49| . 
and in others both shear-thinning and shear-thickening 
are present |T7lil5] . 



The calculation of H17|) is not as direct as it seems. 
First of all it should be noticed that in the fluid state the 
viscosity scales as T^^^ [501150] . The reason behind this Fig. ig. r; normahzed by Tx^^^, friction case. The Newtonian 
dependence is that in these cases, the grain movement is regime is clearer here. The dashed line is explained in the text, 
a consequence of the application of shear, which inputs 




energy in the system. This energy input, for instance, in- 



In flg. E|we plot rj/T^J"^ for frictionless packings. We 



creases the kinetic energy of the grains. Therefore, the see that at low shear rate and low volume fractions we have 
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approximately Newtonian behavior, already expected by 
the result from the shear stress, fig. El This is the effect of 
the high velocities imposed on the system at this regime. 
Particles move so fast that they do not seem to "see" each 
other and then the whole viscosity drops. At higher vol- 
ume fractions we clearly observe a shear-thinning behav- 
ior. We show a fit to the last 5 points for the (f) = 0.7251 
curve. The slope is inside the range reported for many soft 
materials [Bll20n49| . 

These regimes are more clearly observed in the friction 
case, fig. El At lower densities and smaller shear rates, 
we have a Newtonian behavior with a close to zero. The 
situation is a bit different at higher densities and low shear 
rates. The slope of the fit (dashed line) for the first 5 points 
in the <f> = 0.7251 curve is outside the known range for 
complex fluids. Again, friction may be playing a major role 
and more investigation is needed to elucidate this matter. 
It is important to check the behavior in the fluid regime, 
states with = 0.5999, 0.6104 and 0.6318. In the first 
two of them, we have a Newtonian behavior as expected 

1/2 

for r]/T^ . The third one displays shear-thinning since it 
is jammed, see figs. 1151 and 1161 The results for the other 
states should be taken carefully, especially regarding the 
values of the exponents a. But the trend is correct, the 
system should display shear-thinning at high shear rate 
as a result of jamming. 

6 Conclusions 

In this paper we study the effective temperature of sheared, 
bi-disperse, dense granular system in three dimensions as 
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a function of the shear rate and volume fraction. We con- 
sider systems in which grains interact by purely repulsive 
forces, frictionless case, and others that they interact by 
normal and tangential forces, friction case. 

The effective temperature Q is shown to be indepen- 
dent on the shear rate at volume fractions of the order, and 
above, of = 0.68 for the frictionless case, and (p — 0-63 
for the friction case. This effective temperature is shown 
to be equal to both large and small particles, indicating 
that such parameter is at least consistent with a zero-th 
law of thermodynamics for grains. It is shown that T^ff 
has some dependence on the observable as we had only ap- 
proximate agreement with a temperature measured from 
the ISF. However, we present enough evidence pointing to- 
ward the fact that the effective temperature is an intrinsic 
property of the system and, along with experimental re- 
sults 1 13) , that a thermodynamics of granular materials is 
feasible. 

The fact that slow relaxation modes can be character- 
ized by a temperature raises the question of the existence 
of a form of ergodicity for the structural motion, allowing 
a construction of a statistical mechanics ensemble for the 
slow motion of the grains. This argument leads us back 
to the ideas of the thermodynamics of jammed states. In 
parallel to these dynamical measurements, the same infor- 
mation is drawn from the system by a flat statistical av- 
erage over the jammed configurations. Once all the static 
configurations have been visited by the system, the com- 
pactivity Xe can be calculated from the statistics of the 
canonical ensemble of the jammed states. The logarithm 
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of the available configurations at a given energy and vol- 
ume reveals the entropy, from which the compactivity is 
calculated. Our explicit computation shows that the tem- 
perature arising from the Einstein relation ^ can be un- 
derstood in terms of the configurational compactivity Xe 
arising from the statistical ensemble of jammed states. 

We also investigate the role played by the kinetic en- 
ergy in the system. The granular temperature does not 
share the properties of Tgff. It has a power-law depen- 
dence on the shear rate with exponents that increase with 
volume fraction. This is a direct consequence of the en- 
during contacts between grains. 

The jamming transition is observed in the behavior of 
other observable we studied, namely shear stress and ap- 
parent shear viscosity. The states with volume fractions 
higher than <f) = 0.68 and (p ~ 0.63 in the frictionless and 
friction cases develop a yield stress. This is clear indica- 
tion that the system jammed. Moreover, we showed that 
our system obeys the Herschel-Bulkeley law (I16II . We ob- 
served that the exponent is density dependent (in analogy 
with other soft materials, where it is proportional to the 
external bath temperature). This parameter tends, at high 
shear rate and low volume fraction, to 1, in the friction- 
less case, where the system behaves like a power-law fluid 
|19) . We did not observe this in the friction case. In fact at 
lowest volume fraction, the system displays two different 
exponents, and no yield stress. We attribute this to the 
possibility that friction between particle could invalidate 
H16|) for grains in some range of (j). 
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The apparent viscosity, normalized by the square root 
of the granular temperature, displays the predicted behav- 
ior of shear-thinning observed for others complex fluids at 
high volume fractions |19II20II3UII49| . It depends as a power 
law on 7 with an exponent which depends on the volume 
fraction, and approaches —1.0 as (j) increases in the fric- 
tionless case. This exponent is characteristic of complex 
fluids in jammed states. The system shows a Newtonian 
behavior at low volume fraction and shear rate. In the 
friction case, the Newtonian behavior is even clearer, but 
we measured an exponent in eq. H18|) which is outside the 
known range for complex fluids, and this may be, again, 
the presence of friction. However, the trends we obtained 
are important. We observe a shear-thinning behavior for 
the system as long as it stays jammed. All these results 
led to the conclusion that by applying shear to a granular 
packing we can induce a jamming transition (observed for 
the friction case at volume fraction of = 0.6318) or we 
can unjam a system (observed for the frictionless case at 
volume fraction oi cj) ~ 0.6428). 
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